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Estimating  the  frequencies  of  signals  in  a  noisy  environment  has  numerous 
applications  in  digital  signal  processing.  In  December  1980,  Golub  and  Van 
Loan  proposed  a  spectral  estimator  called  the  Total  Least  Squares  (  TLS  ) 
technique  which  is  a  refinement  of  the  Least  Squares  (  LS  )  technique.  In  this 
thesis,  we  first  describe  the  origin  of  the  TLS  technique  and  its  applications  to 
frequency  estimation.  Furthermore,  we  present  a  numerical  implementation  for 
resolving  two  damped  /  undamped  closely-spaced  sinusoidal  signals  in  white 
noise. 

Next,  we  introduce  TLS  extensions  such  as  the  Constrained  Total  Least 
Squares  (  CTLS  )  technique  and  the  Linear  Constraint  Total  Least  Squares 
(  LCTLS  )  technique.  The  CTLS  addresses  the  case  where  the  noise 
components  are  related  and  the  LCTLS  addresses  the  case  where  one  desires 
to  resolve  between  two  narrowband  signals  close  in  frequency,  one  of  which  is 
known.  Finally,  we  present  a  numerical  implementation  of  the  Recursive  Total 
Least  Squares  (  RTLS  )  technique  and  apply  it  to  the  case  of  a  signal  with  a 
fixed  frequency  together  with  a  signal  with  a  time-varying  frequency. 
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I.  INTRODUCTION 


Estimating  the  frequencies  of  signals  in  a  noisy  environment  has  numerous 
applications  in  digital  signal  processing.  Various  techniques  have  been  proposed 
to  solve  that  problem  such  as  the  Pisarenko  algorithm  [1,  p.  623]  and  the 
MUSIC  algorithm  [1,  p.  626].  These  two  algorithms  adopt  the  eigen- 
decomposition  technique  to  obtain  the  solution  from  the  partitioned  signal  and 
noise  subspaces.  In  this  thesis,  we  present  a  subspace  technique  for  frequency 
estimation  called  the  Total  Least  Squares  algorithm.  Two  extensions,  the 
Constrained  Total  Least  Squares  algorithm  and  the  Recursive  Total  Least 
Squares  algorithm,  are  also  introduced. 

In  Chapter  II,  we  first  introduce  the  ordinary  Least  Squares  (  LS  ) 
algorithm  [1,  p.  518].  The  LS  algorithm  is  used  to  solve  the  overdetermined  set 
of  linear  prediction  equations  Ax=b,  where  A  and  b  are  nwi  and  mxl  matrices 
depending  on  data,  and  x  is  an  array  of  parameters.  In  other  words,  the  LS 
algorithm  computes  the  solution  x^  so  that  the  Euclidean  norm  II  Ax-b  II  is 
minimum.  Note  that  the  LS  algorithm  considers  the  errors  or  perturbations  in 
the  observation  vector  b  only.  However,  in  practice,  we  also  need  to  consider 
the  errors  or  perturbations  occurring  in  the  data  matrix  A.  Thus,  we  consider 
the  Total  Least  Squares  (  TLS  )  algorithm  [2]  which  is  an  extension  of  the 
classical  LS  algorithm  and  has  been  applied  to  frequency  estimation  in  signal 
processing  to  locate  signals  in  the  presence  of  additive  noise.  This  algorithm 
can  be  viewed  as  a  generalized  LS  algorithm  and  it  has  better  performance 
than  the  ordinary  LS  algorithm.  However,  it  is  computationally  more  expensive. 
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We  obtain  the  TLS  solution  xTLS  by  simultaneously  minimizing  the 
perturbations  AA  and  Ab  which  exist  in  the  data  matrix  A  and  the  observation 
vector  b,  respectively. 

When  there  exists  some  linear  dependence  algebraic  relationship  among  the 
correlated  noise  perturbation  components  in  A  and  b,  the  TLS  algorithm  has  to 
be  modified  to  take  this  dependence  into  account.  The  modification  leads  to  the 
Constrained  Total  Least  Squares  (  CTLS  )  algorithm  [6].  To  be  more  specific, 
the  CTLS  algorithm  results  from  an  unconstrained  minimization  problem  o'er  a 
small  set  of  variables.  It  uses  the  fact  that  the  elements  in  perturbations  AA 
and  Ab  are  algebraically  correlated.  By  employing  a  new  matrix  F  with  a 
smaller  dimension  and  a  white  noise  random  vector  w,  it  minimizes  the 
Frobenius  norm  of  AC  to  derive  the  solution  xCTLS.  Next,  we  introduce  the 
LCTLS  algorithm  [10].  This  algorithm  can  be  considered  as  an  extension  of 
the  TLS  algorithm  and  a  special  case  of  the  CTLS  algorithm.  It  is  used  to 
resolve  closely-spaced  frequencies  with  a  given  known  frequency  component. 
An  example  and  a  summary  of  the  main  steps  needed  to  be  implemented  the 
LCTLS  are  presented. 

Chapter  III  presents  a  numerical  implementation  of  the  TLS  algorithm  for 
resolving  two  closely-spaced  damped  /  undamped  sinusoids.  We  investigate  the 
effects  of  choosing  different  signal  phase  angles,  signal-to-noise  ratios  (  SNR  ) 
and  prediction  filter  orders.  In  addition,  we  explore  the  estimation  behavior  of 
the  algorithm  due  to  incorrect  singular  subspace  partitions. 

In  Chapter  IV,  we  introduce  the  Recursive  Total  Least  Squares  (  RTLS  ) 
algorithm  [11].  The  algorithm  is  used  in  this  thesis  to  estimate  the  frequencies 
of  two  complex  sinusoidal  signals  in  white  noise.  We  adopt  a  spherical 
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subspace  updating  [12]  technique  and  use  an  eigendecomposition  method  to 
obtain  the  solution  xRTLS.  Four  parameters  are  varied:  the  signal-to-noise 
(  SNR  )  ratio,  the  prediction  filter  order,  the  forgetting  factor  and  the  range  of 
the  frequencies  involved.  We  also  study  the  algorithm  performances  obtained 
by  using  the  signal  or  the  noise  subspace  to  estimate  the  frequencies.  Last, 
Chapter  V  presents  conclusions. 
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II.  THE  TOTAL  LEAST  SQUARES  TECHNIQUE 
AND  ITS  EXTENSIONS 

A.  INTRODUCTION 

Resolving  two  closely  spaced  sinusoidal  signals  in  a  noisy  environment  is  a 
problem  of  importance  in  digital  signal  processing.  This  problem  becomes 
rather  difficult  when  the  number  of  data  samples  is  small  or  the  signal-to-noise 
ratio  (  SNR  )  is  low.  Here,  we  adopt  the  Total  Least  Squares  (  TLS  )  approach 
to  solve  a  Linear  Prediction  Equation  (  LPE  )  to  improve  the  resolution.  Note 
that  the  performance  of  the  algorithm  deteriorates  when  SNR  is  low,  say  below 
3dB. 

The  TLS  technique  is  the  generalization  of  the  ordinary  Least  Squares 
( LS  )  technique  [1]  and  it  is  often  used  to  solve  an  overdetermined  set  of  linear 
equations  obtained  from  noisy  observations.  When  there  are  more  equations 
than  unknowns,  the  system  Ax=b  is  said  to  be  overdetermined.  The  TLS 
method  takes  into  account  the  errors  that  occurs  in  the  data  matrix  A  as  well 
as  in  the  observation  vector  b  and  try  to  attenuate  these  noise  effects  on  both 
sides  simultaneously.  Further  details  are  presented  later  in  this  section. 

B.  LEAST  SQUARES  (  LS  )  ALGORITHM 

I.  Introduction 

First,  recall  the  ordinary  Least  Squares  (  LS  )  technique.  For  the  LS 
algorithm,  the  data  matrix  A  is  assumed  to  be  free  of  error  and  all  errors  are 
confined  to  the  observation  vector  b.  However,  this  assumption  is  not  always 
realistic  since  many  errors  (  perturbations  /  interference  )  like  sampling  errors. 
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modeling  errors  and  instrument  errors,  etc.  may  result  in  inaccuracies  or 
uncertainties  in  A. 

2.  Algorithm  Derivation 

We  obtain  the  LS  solution  x  by  minimizing  the  Euclidean  norm 
II  Ax  -  b  II.  Thus,  it  is  equivalent  to  the  following  expression: 

min|jAb|  (2.1) 

x,Ab 


subject  to 


A  x  =  b  +  Ab  (2.2) 

which  leads  to 

x,  =  (  AH A)1  AHb  =  A+b  (2.3) 

where  A+  is  the  pseudoinverse  of  A,  AHA  is  the  correlation  matrix.  The 
superscript  notation  '  h  '  denotes  the  complex  conjugate  transposition.  Note 
that  the  LS  algorithm  works  well  when  considering  the  noise  only  in  the 
observation  vector  b. 

C.  TOTAL  LEAST  SQUARES  (  TLS  )  ALGORITHM 
1.  Introduction 

First,  we  represent  the  forward  linear  prediction  equation  Ax=b  in  a 
matrix  form: 
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(2.4) 


So 

Si 

•••  Sp-1 

'xP  ■ 

Sp 

Si 

S2 

•••  Sp 

•  • 

•  • 

Xp_i 

= 

Sp+i 

Sn-P-1 

Sn-p 

•  • 

*•*  Sn-2. 

.  X,  . 

_Sn-i_ 

4j  is  the  observed  noisy  signal,  P  is  the  prediction  filter  order  and  is  usually 
larger  than  the  number  M  of  exponentials  but  smaller  than  the  number  N  of  data 
samples. 

If  the  observed  data  is  noisy,  from  a  statistical  viewpoint  the  LS 
solution  will  no  longer  be  optimal  because  it  can  be  biased  and  suffer 
covariance  from  the  accumulation  of  noise  errors  in  AHA.  To  alleviate  this 
problem,  we  adopt  the  Total  Least  Squares  (  TLS  )  technique  [2].  It  has  been 
shown  to  be  equivalent  to  the  minimum  norm  method  [3]  and  it  is  used  to 
diminish  the  bias  accrued  in  AHA  or  to  remove  the  noise  components  by  using 
a  perturbation  analysis  on  A  and  b  of  the  smallest  2-nonn  that  makes  the 
system  equations  consistent.  The  TLS  solution  can  be  obtained  by  minimizing 

II  AC  II  =  II  [  AA|Ab  ]  II  (2.5) 


subject  to 


(  A  +  AA  )  x  =  b  +  Ab 


(2.6) 


which  leads  to 

xtls  =  (  A"  A  •  )'1aH|>  (2.7) 

where  AA  and  Ab  are  two  independent  noise  perturbations  and  a2  is  the 
minimum  eigenvalue  of  [  Alb  ]H[  Alb  j. 
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Note  that  we  use  a  zero  mean  white  identically  independent 
distributed  (  i.i.d. )  Gaussian  noise  in  our  analysis  of  the  TLS  algorithm.  When 
the  noise  is  not  white  but  colored  and  still  uncorrelated  with  the  signals,  a 
whitening  filter  ought  to  be  applied.  This  kind  of  filter  is  capable  of  transforming 
a  stationary  discrete-time  non-white  process  at  the  input  into  a  white 
uncorrelated  data  sequence  at  the  output.  The  whitening  procedure  takes  the 
non-white  noise  vector  w  and  computes  R=E[ww*]«ZZ*,  using  a  Cholesky 
decomposition.  This  leads  to  the  definition  of  a  new  white  noise  vector  5  where 

5  =  Z  *w.  (2.8) 

The  whitening  procedure  can  be  performed  by  using  a  prediction-error  filter 
which  offers  a  linear  prediction  equation  same  as  Ax=b.  Basically,  the 
prediction  is  dependent  on  the  presence  of  correlation  between  adjacent 
samples  of  the  input  non-white  stochastic  process.  Then,  we  can  successively 
decrease  the  correlations  by  increasing  the  prediction-error  filter  order  high 
enough  until  it  ultimately  gets  to  a  point  at  which  the  output  sequence  is 
composed  of  (  white  )  uncorrelated  samples  [4,  p.  216]. 

The  TLS  solution  xTLS  is  obtained  by  using  information  obtained  from 
the  singular  value  decomposition  (  SVD  )  of  the  matrix  C=[  Alb  ]  with  a 
dimension  (N-P)x(P+l),  i.e., 

C  =  UZVH  (2.9) 

where  U  and  V  are  unitary  square  matrices  with  a  dimension  (N-P)x(N-P)  and 
(P+l)x(P+l),  respectively;  Z  is  a  rectangular  diagonal  matrix  with  a  dimension 
(N-P)x(P+l)  consisting  of  the  singular  values  in  decreasing  order  with  the  largest 
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one  on  the  upper-left-handed  comer.  The  matrices  U,  V  and  £  can  be 
partitioned  as 


u  =  [  m  i  u2 1  •  •  •  1  %-p  ] 

(2.10) 

V  as  [  Vj  1  V2  1  •  •  •  1  Vp+1  ] 

(2.11) 

and 

£  =  diag  ( <Tj,  •  •  • ,  Gp^  Op+j  ) . 

(2.12) 

The  columns  of  the  matrices  U  and  V  are  the  unitary  eigenvectors  of  AAH 
and  AhA  respectively. 

2.  Algorithm  Derivation 

The  TLS  algorithm  described  above  can  be  summarized  into  the 
following  equations  [5]  : 

( [  Alb  ]  +  [  AAlAb  ] )  =  0  (2.13) 

or 

(  C  +  AC  )  z  =  0  (2.14) 

where 

C  =  [  Alb  ],  AC  =  [  AAlAb],  z  =  (2.15) 

From  the  homogeneous  equation  (2.13),  we  see  that  the  TLS  problem 
finds  the  minimum  norm  of  the  perturbation  matrix  AC  such  that  (  C+AC  )  is 
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rank  deficient.  In  other  words,  the  TLS  solution  can  be  formulated  as  seeking  a 
solution  vector  x  [6]  to  the  following  problem 


min||AC|£ 

AC,X  * 


(2.16) 


subject  to 


(  C  +  AC  ) 


=  0 


(2.17) 


where 

||ACg  =  tr{AC*AC},  (2.18) 

II  •  llF  denotes  the  Frobenius  norm  and  "  tr  "  is  the  t race  of  a  matrix. 

Next,  the  matrix  C  is  decomposed  in  the  following  form 


C  =  [Uj 


1 — 

0 

rv,Hi 

- 1 

CM 

w 

— 1 

=  <s 
> 

_ 1 

(2.19) 


where  Ui,  £i  and  Vi  correspond  to  the  signal  subspace;  U2,  £2  and  V2 
correspond  to  the  noise  subspace.  We  separate  the  signal  subspace  from  the 
noise  subspace  by  using  the  number  of  sinusoids.  Two  different  cases  need  to 
be  considered.  The  first  one  considers  the  case  where  the  singular  values  are 
all  distinct.  The  second  one  considers  the  case  where  some  of  the  singular 
values  are  multiple.  For  the  first  case,  we  usually  assume  that  the  last  element 
ap+1  in  the  diagonal  matrix  £  is  the  smallest  singular  value  and  it  corresponds  to 
the  last  vector  vp+1  in  the  unitary  matrix  V.  Then  we  get 
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(2.20) 


V(P+D, 

V(P+D, 

where  (vP+1)p+1  is  the  last  element  of  the  last  singular  vector  vP+1  and  is 
assumed  to  be  non-zero  for  normalization.  If  (vP+1)p+1  is  equal  to  zero,  a  TLS 
solution  does  not  exist.  For  the  second  case,  we  assume  that  the  P-M+l 
smallest  singular  values  are  equal,  i.e., 

°i  -°2  >Om+i  =**,=<yp+if  (2.21) 

and  set  V2  to  be  a  column  partition  of  V  where 

V2  -  (  VM+1,  •  •  • ,  vp+1 )  (2.22) 


or 


(2.23) 


where  cT  represents  the  first  row  of  the  matrix  V2. 

Next,  we  introduce  the  Householder  (  reflection  )  matrix  Q  which  is 
used  to  solve  the  TLS  problem.  It  has  the  property  to  zero  out  specific  entries 
in  a  matrix  or  vector.  Usually,  it  is  defined  as  [4,  p.  428]: 


w 

Q  =  Ip  -2— v— 

V  v  v 


(2.24) 
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where 


I  is  an  mxm  identity  matrix, 

v  is  an  mxl  complex  vector  which  corresponds  to  the  complex  vector  Vjin 
V2  defined  in  (2.22)  and  has  the  Euclidean  norm 

l 

|vfl  =  (vHv)2i  (2.25) 

The  vector  v  can  be  defined  below  for  analysis: 

v  =  c*  -  a*eP  =>  vH  =  cT  -  aTeJ  (2.26) 

a‘=±A|cI  =>  a  =  ±r£r|cl  (2.27) 

|Cp I  |Cp| 

ej  =  [0,  — ,  0,  l]lxP=>  Vp  =cP-a  (2.28) 

For  any  given  vector,  we  can  construct  a  Householder  matrix  so  that  all  the 
energy  is  compressed  into  a  selected  partitioned  column  vector  [7].  This  leads 
to  the  following  TLS  solution: 


V2Q  = 


'0 

w 


a 

y 


(2.29) 


cTQ  =  [  0,  •  •  • ,  0,  a  ]. 


(2.30) 


Now  we  express  V2Q  to  extract  the  last  column,  which  is  used  in  the 

solution: 
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V2Q  =  V2  [  q1 1  q2  i  q3 1  •  ■  •  I  qp  ] 


(2.31) 


and  let 


V2qp  = 


(2.32) 


Note  that 


vHv  =  (cT  -  aTe?  )(c*  -  a*eP ) 

=  |c||2  -  ct*cP  -  acp  +  ||a|p  (va*cP  =ac*P) 


=  |cl2-2a*cp+||af 


=  ||cf -2^j-||cO+^l-IcI2  (v|a|2  =  aa‘) 


-2  ||af-^f|cl  (v||o!2=Ic|2=cHc) 

^  IcpI  J 

=  2(aa*-a‘cP). 


(2.33) 


The  last  column  of  the  matrix  Q  is  given  by 


2v  v 
<ip  =  ep-Ti£r 


2(cP-a)  /  •  *  \ 

=  -a'e) 


*  * 


*  * 


Cp  a*(a-cP) 

c*(a-cP)  _  c* 
a*(a-cP)  a*. 


cPc  -cPa  ep-ac  +aa  eP 


(2.34) 
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Thus, 


v'  V 
y=v2qp=^- 


(235) 


y_  v2c*_  v2c* 


XTLS  “ - —  " 


a  aa 


(2.36) 


Equation  (2.36)  corresponds  to  the  noise  subspace  case  since  V2  and  c  are 
obtained  from  V2  which  is  described  in  (2.19).  Similarly,  we  can  derive  the 
solution  vector  from  the  signal  subspace.  First,  we  use  the  identity 


r„H  ^tnf  .7*h  "j 

8  c  8  vi 

_v;  v2J[c  v2hJ~ 


(2.37) 


which  leads  to 


gHg  +  cHc  =  l 
v;g+v2c=o 


cHc  =  l-gHg 

v2c=-v;g  . 


(2.38) 


Next,  we  substitute  (2.38)  into  (2.36)  and  obtain  the  alternative  TLS  solution 
corresponding  to  the  signal  subspace  case: 


x  V'g 

XTLS  -  1  -  gHg 


(2.39) 


where  Vj  and  g  are  obtained  from  Vi  which  is  also  described  in  (2.19). 
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D.  CONSTRAINED  TOTAL  LEAST  SQUARES  (CTLS) 

ALGORITHM 

1.  Introduction 

When  the  noise  components  contained  in  the  augmented  matrix  C 
are  algebraically  correlated  or  linearly  dependent  (  so  called  constraints  ),  the 
TLS  technique  is  no  longer  an  optimal  frequency  spectral  estimator.  In  this 
case,  we  have  to  reformulate  the  problem  to  reduce  the  dimensionality  of  the 
disturbances.  Thus,  the  TLS  needs  to  be  reformulated  to  reduce  the 
dimensionality  of  the  noise  entries.  Thus,  the  Constrained  Total  Least  Squares 
(  CTLS  )  algorithm,  proposed  by  Abatzoglou  et  al  [6],  can  be  viewed  as  the 
natural  extension  of  the  TLS  algorithm.  In  addition,  the  CTLS  can  be  shown  to 
be  equivalent  to  the  Maximum  Likelihood  (  ML  )  algorithm.  The  CTLS  solution 
can  be  obtained  from  a  constrained  state- space-parameter  ML  estimator 
together  with  a.  Newton  algorithm  [8]. 

2.  Algorithm  Derivation 

In  order  to  remove  the  correlated  noise  perturbation  components 
among  A  and  b,  we  derive  the  CTLS  algorithm  from  an  unconstrained 
minimization  problem  over  a  small  set  of  variables.  Let  us  redefine  AA  and  Ab 
both  as  linearly  constrained. 


A  A  =  [  FiW  1  F2w  1  •  •  •  1  Fpw  ] 

(2.40) 

Ab  =  FP+1w 

(2.41) 

Then  the  augmented  error  matrix  AC  is  formed  by 

AC  =  [  AAlAb  ]  =  [  Fjw  I  F2w  I  •  •  •  I  FPw  |  FP+1w  ]  (2.42) 
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or 


AC  =  [  ACj  I  •  •  •  I  ACP  |  ACP+1  ]  =  Fw,  (2.43) 

F  =  [  Fj  I F2  I  •  •  •  I  FP  |  FP+1  ]  (2.44) 

where  F  is  a  matrix  whose  elements  are  related  to  those  of  independent 
variables  AC  and  w  is  a  zero  mean  white  noise  random  vector  of  minimal 
possible  dimensionality.  The  reason  for  w  holding  the  minimum  dimension  is  to 
increase  the  number  of  algebraic  relationships  the  entries  of  AC  can  satisfy  and 
to  accomplish  the  minimization  faster.  So  the  CTLS  solution  can  be  formulated 
as  follows: 


min|w||2  (2.45) 

w,x 


subject  to 


( C  +  AC  ) 


(2.46) 


Next,  we  consider  the  quadratic  constraint  equation  in  (2.46)  which  is  expanded 
into: 


(2.47) 
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where 


AC  j  =  [  FjW  I  F2w  I  •  •  •  I  Fpw  |  Fp+1w  ] 


Let 


=  ZxiFiwFp+iw 

i=l 


=^ZxiFi-Fp+ijw 


then  we  have 


H*s  5>iFrFP+i> 

i=l 


0+h*w=o  =*  w=_h«c(xi) 


where  Hx  is  the  pseudoinverse  of  Hx  defined  as 


Hj*  h;(h„h;)‘. 


For  any  w,  the  solution  vector  x  can  satisfy  (2.45)  and  (2.50);  i.e., 


min  w 


=(xiJc-Hi*Hic(x1). 


Next,  let  us  define  the  analytic  function  J(x) 


«*)■(*)  c'(h;)‘h;c[x] 


(2.48) 

(2.49) 

(2.50) 

(2.51) 
we  have 

(2.52) 

(2.53) 
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X  =  (  Xj,  .  .  .  ,  Xp  ). 


J(x)  can  be  simplified  to  the  following  expressions  by  substituting  (2.51)  into 
(2.53). 

c'hi(h,h;)'1h;(h,h;)1cP 
= P  c*[(h,h;)1h,h;](hih;)1c^^ 

=  '*)  (254) 

Hence,  we  conclude  that  the  CTLS  solution  is  obtained  by  minimizing 
J(x).  From  (2.45)  and  (2.52),  we  try  to  find  out  the  minimum  of  J(x)  with  a 
respect  to  x.  It  is  necessary  to  compute  the  first-order  complex  partial 
derivative  or  complex  gradient  and  set  the  derivative  to  be  zero;  namely. 


aj 

3x 


(2.55) 


Due  to  the  complex  data,  Abatzoglou  et  al  [6]  defined  a  complex 
version  of  the  Newton  method  which  uses  the  second-order  complex  partial 
derivative  of  J(x)  to  update  the  values  of  x.  The  Newton  iteration  (  recursion  ) 
is  briefly  described  by 

x(n  + 1)  =  x(n)  +  (AB-1A  -  b)"1  (a  -  AB_1a)  (2.56) 
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where 


a  =  =  f  •  •  •- — 1  =  complex  gradient  of  J 

3x  3xp  J 

1  (  92J  32jt  ^ 

A  =  — I  — j  +  -  j-  =  unconjugated  complex  Hessian  of  J 

2  ^  ox *  ox*  J 

_  a2j 

»  =  ^-0'"  =  conjugated  complex  Hessian  of  J.  (2.57) 

Hence,  by  referring  to  the  derivation  in  [6],  we  obtain  the  following 
closed-form  solutions  for  a,  A  and  B: 


a  =  (h  *  B)t 

A  =  -F*  Hi(H,H*  r1  B-(F*  H^HX)'1  B)t 

b = [b*  (hxh;  r1  b]t+ f*  [h;  (hxh;  r1  hx  -  n  f  (2.58) 


where 


h=(HxH;r1c[^_iJ 

B  =  CIP+1>P  -[FjHjhl—lFpHjh] 

G  =  [F1h|”*|Fph].  (2.59) 


The  Newton  method  is  implemented  by  starting  with  an  Iterative 
Quadratic  Maximum  Likelihood  ( IQML  )  algorithm  [9].  This  algorithm  solves 
a  quadratic  constraint  minimization  problem  and  derives  the  solution  vector 
from  the  coefficients  of  the  linear  prediction  polynomial  which  correspond  to  the 
parameters  of  the  signals  in  a  fast  convergent  way.  We  summarize  its  main 
steps  below: 
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(d) .set  n=n+l  for  real-time  iteration 

(e) .check  the  convergence:  II  x(n-l)  -  x(n)  II  <  e 

(f) .identify  the  roots  of  the  prediction  vector  polynomial  x  which  represent  the 
parameters  of  the  received  signal. 


19 


E.  LINEARLY  CONSTRAINT  TOTAL  LEAST  SQUARES 
(LCTLS)  ALGORITHM 

1.  Introduction 

In  this  section  we  introduce  a  special  case  of  the  CTLS  algorithm 
called  the  Linearly  Constrained  Total  Least  Squares  (  LCTLS  )  technique.  The 
LCTLS  technique  is  used  to  locate  and  resolve  closely  spaced  frequencies  in 
the  presence  of  a  given  known  frequency  component  [10]. 

Recall  that  the  TLS  method  solves  the  unconstrained  problem  : 

(  A  +  AA  )  x  =  b  +  Ab.  (2.65) 

Using  (2.13)  to  (2.18),  we  compute  the  minimum  norm  of  the  total  perturbation 
matrix  AC  which  makes  (  C  +  AC  )  rank  deficient  and  generates  a  vector  z 
defined  from  a  null  space.  Equation  (2.65)  has  a  solution  provided 

(  b  +  Ab  )  e  Range  (  A  +  AA  ).  (2.66) 

In  this  case,  Zn+i*0  and  we  can  solve  for  x.  If  z„+i=0,  the  TLS  problem  has  no 
solution.  In  addition,  the  non-zero  last  element  Zn+i  is  used  to  normalize  the 
vector  z;  i.e  .  we  have 

•(£)■(-■)  - 

2.  Algorithm  Derivation 

To  solve  the  LCTLS  problem,  we  first  look  at  the  constraint  equation: 

¥  x  =  T  (2.68) 
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where  'F,  x  and  T  have  dimensions  px(n+l),  (n+l)xl  and  pxl,  respectively. 
Equation  (2.69)  can  also  be  expressed  as  below: 

[  'P  r  ]  j  =  [  ¥  T  ]  z  =  0.  (2.69) 

Let  us  introduce  a  physical  example  to  define  the  matrix  [  *¥  T  ].  First,  we 
consider  a  signal  £(n)  composed  of  two  complex  sinusoidal  signals  in  white 
noise: 

£(n)  =  A1eino>1  +  A2e-’ne0j  +w(n),  n=  1,  •••  ,N  (2.70) 

where  £(n)  is  the  observed  noisy  data;  w(n)  is  the  white  Gaussian  noise;  Ai  is 
the  amplitude;  ©i  is  the  angular  frequency  and  N  is  the  data  length. 

Assume  that  0)i  is  known  and  0)2  is  unknown.  Create  a  signal 
source  vector  X(©)  defined  as: 

JtH(©)  =  [1,  ej®,  ej2<0,  ej3a> . ejn(0]>  (2.71) 

We  make  A.(oo)  satisfy  the  constraint  equation  in  (2.69),  which  leads  to 


=  0. 


(2.72) 


The  constraint  forces  e*0'  to  be  one  of  the  roots  of  the  prediction-error 
polynomial  e(o))  defined  as: 


e(©)  =  XH(co) 


-1 


=  0. 


(2.73) 
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Then,  the  angle  of  the  next  closest  root  on  the  unit  circle  is  chosen  as  the 
estimate  of  the  unknown  angular  frequency  u>2. 

Next,  considering  the  TLS  equation  given  in  (2.65)  and  constraint 
given  in  (2.69)  leads  to 

(a-TaA  b+rAb)(-l)  =  0-  <2-74> 

Next,  using  a  SVD  partition  method  to  obtain  the  solution.  To  ensure  the 
constraint  equation  in  (2.75)  to  be  satisfied,  we  have  to  constrain  the  vector  z  to 
fall  in  the  null  space  of  [  T  ]  which  can  be  further  factored  out  using  QR 
decomposition  into 

[  r  ]  =  [  Q,  0  ]AH  =  QjA?  (2.75) 

-tD  (2-76) 

where  Qi  is  an  upper  triangular  matrix  with  a  dimension  (pxp)  and  AH  is  a 
unitary  matrix  with  a  dimension  (n+l)x(n+l).  By  comparing  the  matrix 
dimensions  from  (2.69)  and  (2.75),  we  see  that  z  should  be  the  product  of  the 
orthonormal  basis  A2with  a  dimension  (n+l)x(n+l-p)  for  the  null  space  of 
[  r  ]  and  an  arbitrary  non-zero  vector  r\  with  a  dimension  (n+l-p)xl,  i.e.. 


z  =  A2tv 


(2.77) 
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Next,  define  the  matrix  C  as 


C  =  CA2a5 


(2.78) 


with 


C  =  [A!  b].  (2.79) 

The  SVD  of  the  matrix  C  is  performed  and  C  is  partitioned  into  its  signal  and 
noise  subspaces,  which  leads  to: 


_  _  „  «. H  n+l-p  _  «.  *  H 

c=uzv  =  X  °iuivi  (2.80) 


0 
0 

with 

Vi  :  r-dimensional  unconstrained  part  of  the  signal  subspace 

V2  :  s-dimensional  noise  subspace 

V3  :  p-dimensional  constrained  part 
*  « 

£i  =diag(ai,---,ar) 

S2  =on+i-pIs,  s=n+i-p-r 
£3  =0. 


rH 


V" 


(2.81) 


i=l 


or 


C-jlJj  U2  UjJ 


Ej  0 

0  i2 
0  0 
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Using  the  fact  that 


AjA^  +  A2A2  =  I, 


(2.82) 


we  get 


c=ca1a?+ca2a? 


(2.83) 


Dowling  et  al  [10]  have  suggested  to  enforce  sphericity  in  the 
concerned  partitioned  subspace  to  improve  the  frequency  spectral  estimation 
performance.  The  term  "  sphericity  "  means  that  the  singular  values  in  the 
signal  and  noise  spaces  are  replaced  with  their  mean  values.  For  example, 
assume  that  there  exists  (  s=n+i-p-r  )  equal  singular  values  in  the  noise 
subspace  V2;  i.e.,  tlr  noise  singular  values  are  replaced  with  the  variance  of 
the  random  noise.  By  partitioning  the  truncated  unitary  matrix  V  below  and 
doing  the  same  operations  as  in  (2.37)  and  (2.38) 

v  =  [v,  v2  V,]  =  U  J  jl  (2.84) 


we  obtain  the  minimum  norm  solution  corresponding  to  the  noise  subspace: 


LCTLS  ~  CHC 


(2.85) 
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Further  using  the  following  identity 


which  leads  to  the  following  two  equations: 

gHg+cHc  +  fHf  =  l  =>  cHc  =  l-gHg-fHf 
V1'g  +  V2c  +  V3f  =  0  =>  V2c  =  -Vlg-V3f 


(2.86) 


(2.87) 


we  obtain  the  minimum  norm  solution  corresponding  to  the  signal  subspace: 


XLCTLS  =  i-gHg-f^f 


(2.88) 


Note  that  the  minimum  LCTLS  solution  should  be  based  on  the  subspace  V2  or 
[%  V3 whichever  one  has  the  smaller  dimension. 
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III.  TOTAL  LEAST  SQUARES  TECHNIQUE 
IMPLEMENTATION 


A.  INTRODUCTION 

First,  we  study  the  performance  of  the  TLS  algorithm  to  estimate  the 
frequency  of  sinusoidal  signals  in  noise.  We  assume  the  received  signal  x(n) 
consists  of  M  uniformly  spaced  signal  data  samples  corrupted  by  a  complex 
white  Gaussian  noise  w(n).  The  received  signal  is  expressed  by  the  following 
equation  and  is  observed  over  a  duration  of  N  seconds,  assuming  the  sampling 
frequency  fs=lHz  : 

M 

x(n)  =  Xakej2rtf‘n  +  w(n),  n  =  0,  1,-,  N-l  (3.1) 

k=l 


where 

{  ak  }  :  a  set  of  signal  amplitudes 
{  fk  )  :  a  set  of  true  signal  frequencies 
{  x(n)  }  :  measured  complex-valued  data  samples 
{  w(n)  }  :  random  complex-valued  noise  samples  and  uncorrelated 
with  signals 

Real  and  imaginary  part  of  w(n)  are  assumed  to  have  a  variance  a2. 

We  are  interested  in  estimating  the  signal  frequency  parameters  fk;  the 
parameters  ak  are  chosen  to  have  unit  magnitude  for  the  sake  of  simplicity. 
Note  that  we  use  the  singular  value  decomposition  (  SVD  )  technique  to  solve 
the  linear  prediction  equation  Ax=b  which  matrix  form  is  given  in  (2.4).  The 
TLS  implementation  is  expected  to  decrease  the  noise  effects  occurring  both  in 
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the  data  matrix  A  and  in  the  observation  vector  b  simultaneously.  First  we 
choose  the  order  P  of  the  prediction-error  filter  and  solve  a  P-degree  prediction- 
error  filter  polynomial.  Next,  we  determine  the  signal  parameters  from  the  M 
zeros  which  are  the  closest  to  the  unit  circle.  The  other  P-M  zeros  are  called 
extraneous  zeros;  they  are  approximately  uniformly  distributed  inside  the  unit 
circle.  Note  that  the  choice  of  optimal  prediction-error  filter  order  is  important 
since  if  the  filter  order  is  chosen  too  low,  the  closely-spaced  sinusoids  can  not 
be  resolved.  In  addition,  if  the  filter  order  is  chosen  too  high,  the  noise  zeros 
may  be  mistaken  for  the  true  signal  zeros  as  the  noise  zeros  become  very  close 
to  the  unit  circle.  Thus,  we  have  to  select  the  "  optimal  "  prediction  filter  order 
carefully. 

B.  ALGORITHM  IMPLEMENTATION 

The  parameters  used  in  the  simulation  are  <j>i=— ,  <J>2=0;  ^  and  f2  represent 

4 

the  frequency  locations;  w(n)  are  independent  identically  distributed  (  i.i.d.  ) 
white  complex  Gaussian  noise  samples  with  zero  mean  and  variance  a2  for 
both  real  and  imaginary  parts.  The  signal-to-noise  ratio  (  SNR  )  of  the  received 
noisy  signal  is  given  by  the  peak  amplitude  of  sinusoids  to  noise  level  and 

defined  as  10*log(— )dB.  The  sampling  frequency  is  set  to  be  1Hz.  The 
2a 

frequency  separation  is  chosen  smaller  than  the  reciprocal  of  the  observation 
time  in  order  to  resolve  these  two  close-together  frequencies  and  raise  the 
resolution.  Simulations  are  performed  using  an  ensemble  of  30  independent 
trials.  For  each  trial,  a  data  block  of  50  (  N=50  )  data  samples  is  used  to 
resolve  two  (  M=2  )  sinusoids.  The  frequency  estimates  are  computed  by 
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choosing  the  two  roots  of  the  prediction-error  filter  transfer  function  or  the  two 
zeros  of  the  prediction-error  filter  polynomial  which  are  the  closest  to  the  unit 
circle  in  the  complex  z-plane.  Different  values  of  the  filter  order  P  in  the  range 
of  2~49  can  be  used  to  observe  the  variations  of  estimation. 

C.  EXPERIMENTAL  OBSERVATIONS 
1.  Performance  of  the  Algorithm 

We  investigate  the  performance  of  the  algorithm  with  a  statistical 
quantity  called  the  effective  standard  deviation  (  ESD  )  defined  below  : 


(3.2) 


where 

T  :  the  total  number  of  independent  trials 

A 

f,  :  the  estimated  frequency  at  the  i*  realization 
f  :  the  true  signal  frequency. 

The  ESD  gives  a  measure  of  the  distance  between  the  true  and  the  estimated 
frequency  locations.  In  this  simulation,  we  first  use  fj=0.23Hz  while  f2  =0.27Hz. 

K 

The  data  including  <)>,=— ,  <J>2=0,  P=8,  N=50  and  30  independent  trials  are  used 

4 

to  compute  the  ESD.  The  selection  of  the  two  signal  zeros  is  made  by  creating 
two  sector  regions  designated  from  the  origin  to  point  z1=0.3+j*0.7  and  to  point 
z2=  -0.3+j*0.7,  respectively.  Varying  the  SNR  from  20dB  to  OdB,  we  observe 
that  the  sector  region  for  fi=0.23Hz  at  SNR=0dB  needs  to  be  expanded  by 
changing  Z!=0.4+j*0.7  and  the  sector  region  for  f2=0.27Hz  at  SNR=ldB  needs 
to  be  expanded  by  changing  z2=-0.4+j*0.8.  Here  we  display  six  different  SNR 


28 


cases  which  are  20dB,  12dB,  5dB,  3dB,  ldB  and  OdB.  We  also  observe  that  the 
signal  and  noise  zeros  resolution  becomes  worse  as  the  SNR  becomes  smaller, 
as  shown  in  Figures  3.1.1-3.1.6.  Figure  3.1.6  shows  the  worst  resolution  as  a 
couple  of  noise  zeros  are  on  the  unit  circle  and  could  potentially  be  considered 
as  the  signal  zeros.  In  addition.  Table  1  shows  that  as  the  SNR  becomes 
smaller,  the  estimated  deviations  for  both  signal  frequencies  are  getting  larger. 
The  above  simulation  results  compare  with  those  obtained  in  [5]. 

2.  Effects  of  Signal  Phase  Changes 

Two  sets  of  phases  are  chosen  to  investigate  the  effects  of  phase 

changes:  (1)  <j>i=— ,  <j>2=0,  (2)  <j)i=— ,  <J>2=-— .  The  filter  order  is  P=8;  the 
4  4  4 

number  of  data  samples  points  is  N=50;  30  independent  trials  are  used;  the 
signal  frequencies  are  chosen  as  fi=0.25Hz,  f2=0.27Hz;  the  SNR  values 
considered  are  3dB  and  12dB.  Figures  3.2.1~3.2.4  show  that  phase  changes  do 
not  affect  the  resolution. 

3.  Effects  of  Filter  Order  Changes 

Next,  we  consider  the  effect  of  changes  in  the  prediction-error  filter 
order.  Recall  that  it  is  important  to  choose  the  smallest  value  of  P  which 
achieves  a  good  frequency  estimation  performance.  For  this  experiment,  we 
choose  fi=0.25Hz,  f2=0.27Hz,  SNR=12dB,  N=50,  and  10  independent  trials  are 
used.  Figures  3.3.1  to  3.3.6  show  the  effects  of  changing  the  filter  order.  When 
P=2,  we  have  very  poor  resolution,  as  shown  in  Figure  3.3.1.  Performance 
slightly  improves  when  P=4,  as  shown  in  Figure  3.3.2.  Note  that  increasing  the 
filter  order  P  improves  the  resolution  of  the  estimated  frequencies.  For 
example,  when  choosing  P=8  and  P=18,  we  observe  that  two  sinusoidal  signals 
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are  located  exactly  at  their  desired  positions  on  the  unit  circle  and  the 
extraneous  zeros  are  uniformly  distributed  inside  the  unit  circle,  as  shown  in 
Figures  3.3.3  &  3.3.4.  However,  if  P  is  chosen  too  large,  P=28  or  P=38  for 
example,  we  observe  some  of  the  P-M  extraneous  zeros  of  the  prediction 
polynomial  are  getting  very  close  to  the  unit  circle  which  make  the  distinction 
between  signal  and  noise  zeros  difficult  to  make,  as  shown  in  Figures  3.3.5  & 
3.3.6. 

4.  Effects  of  Errors  in  the  Subspace  Partition  Estimate 

It  is  possible  to  express  the  TLS  solution  in  the  sense  of  minimum 
norm.  By  applying  the  SVD  technique  to  the  linear  prediction  equation,  we  can 
partition  the  unitary  matrix  V  (  for  handling  complex  data  )  into  two  subspaces: 
the  signal  subspace  and  the  noise  subspace.  If  the  partition  is  performed 
correctly,  we  can  obtain  a  good  zeros  resolution  by  using  a  moderate  filter 
order  P  and  high  enough  SNR,  choosing  P=8  and  SNR=12dB  for  example,  as 
shown  in  Figure  3.3.3.  However,  if  the  partition  is  not  performed  correctly;  i.e., 
the  numbers  of  column  vector  of  the  signal  subspace  and  the  noise  subspace 
are  not  estimated  properly,  we  observe  that  the  extraneous  zeros  of  the 
prediction  polynomial  are  randomly  distributed.  The  more  the  signal  and  the 
noise  subspaces  are  incorrectly  chosen,  i.e.,  the  worse  the  partition  is,  the 
worse  the  resolution  to  identify  the  true  sinusoidal  signal  locations  becomes,  as 
shown  in  Figures  3.4.1  &  3.4.2.  Figure  3.4.2  displays  the  worst  degradation  of 
the  performance  since  these  two  subspaces  are  mixed  up  by  6  column  vectors 
under  the  condition  of  SNR=12dB,  P=8,  N=50  and  10  independent  trials. 


30 


5.  Application  to  Damped  Signals 

In  the  following,  we  consider  the  case  of  a  sum  of  two  (  M=2  ) 
damped  exponential  sinusoids  described  as  follows  : 

x(n)=  +  e,-“*"i<2^n+«>)4w(n),  (3.3) 

n  =  0, 1,  •  •  • ,  N-l 

where  angular  frequencies  are  f!=0.25Hz,  f2=0.27Hz;  the  signal  phases  are 
fixed  at  <j>i=— ,  <J>2=0;  the  two  damping  factors  are  chosen  aj=0.8,  a2=0.9;  the 

number  of  data  sample  is  N=50;  the  filter  order  is  P=8;  30  independent  trials  are 
used  for  the  simulation,  and  w(n)  is  still  defined  as  the  complex  white  Gaussian 
noise.  We  observe  that  it  is  quite  difficult  to  identify  the  exact  locations  of  these 
two  damped  sinusoidal  signals  at  low  SNR,  for  example  5dB,  or  at  high  SNR, 
for  example  20dB,  as  shown  in  Figures  3.5.1  &  3.5.2.  Even  if  we  increase  the 
filter  order  up  to  P=18,  the  resolution  has  not  been  improved  still,  as  shown  in 
Figures  3.5.3  &  3.5.4. 
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IV.  RECURSIVE  TLS  TECHNIQUE 
IMPLEMENTATION 


A.  THEORETICAL  BACKGROUND 

1.  Introduction 

As  discussed  before,  the  Total  Least  Squares  (  TLS  )  algorithm  has 
been  shown  to  be  more  accurate  than  the  Least  Squares  (  LS  )  algorithm  when 
the  errors  or  noise  perturbations  exist  on  both  sides  of  a  linear  prediction 
equation  Ax=b.  We  now  introduce  an  adaptive  algorithm  called  the  Recursive 
Total  Least  Squares  (  RTLS  )  algorithm.  In  this  algorithm,  we  perform  the 
updates  by  tracking  the  smaller  one  of  the  signal  subspace  or  the  noise 
subspace  partitioned  from  a  unitary  matrix  containing  complex  data.  It  is 
shown  that  the  RTLS  algorithm  works  better  than  the  Recursive  Least  Squares 
(  RLS  )  algorithm  [4,  p.  480]  and  the  Least  Mean  Squares  (  LMS  )  algorithm 
[4,  p.  302]  in  terms  of  tracking  property  as  well  as  steady-state  tap-weight  error 
norms  [11].  Furthermore,  the  updates  can  also  be  done  using  a  computationally 
efficient  non-iterative  subspace  updating  technique  [12]. 

2.  Algorithm  Derivation 

First,  let  us  give  a  brief  review  of  an  adaptive  process.  An  adaptive 
process  involves  the  automatical  adjustment  of  a  set  of  tap  weights  of  a 
prediction-error  filter;  i.e.,  it  obtains  an  estimate  of  a  desired  response  which 
results  from  the  inner  product  of  a  set  of  tap  input  and  its  corresponding  set  of 
tap  weights  produced  by  the  adaptive  filter.  Then,  it  compares  this  obtained 
estimate  with  the  actual  values  of  the  desired  response  and  generates  an 
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estimation  error.  The  estimation  error  is  fed  back  to  the  adaptive  process  and 
stimulates  the  prediction  filter  in  a  recursive  way.  Therefore,  this  is  a  closed- 
loop  operation  and  the  final  goal  of  the  adaptive  process  is  to  make  the  output 
error  reach  its  minimum  value.  There  are  two  good  examples:  the  LMS 
algorithm  and  the  RLS  algorithm. 

In  the  TLS  algorithm,  we  apply  the  singular  value  decomposition  (  SVD  ) 
technique  to  the  augmented  matrix  spanned  by  the  data  matrix  A  and  the 
observation  vector  b.  However,  in  the  RTLS  algorithm,  we  apply  an  eigen- 
decomposition  technique  to  a  specific  augmented  matrix  C(n)  defined  below  to 
get  the  desired  subspaces  and  determine  the  effective  rank  of  C(n)  from  its 
singular  values.  In  addition,  we  adopt  a  parameter  a,  called  the  exponential 
weighting  factor ,  which  is  similar  to  X  (forgetting  factor  )  in  the  RLS  algorithm. 

C(n)  =  D<”>  [  A  I  b  ]  (4.1) 

where 

D<an)  =  VdUg(an . a0).  (4.2) 

It  is  worthy  to  investigate  the  impact  of  the  forgetting  weight  factor  a  on 
the  tracking  capability.  The  factor  a  is  used  to  ensure  that  the  data  samples  in 
the  distant  past  are  "  forgotten  "  so  as  to  afford  the  possibility  of  following  the 
statistical  variations  of  the  observable  data  when  the  filter  operates  in  a  non- 
stationary  environment.  It  is  designated  to  be  a  positive  constant  close  to,  but 
smaller  than  one.  Further  speaking,  the  reciprocal  of  1-a  can  be  considered  as 
a  measure  of  "  memory  ".  When  a  is  equal  to  1 ,  which  corresponds  to  infinite 
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memory  or  most  unforgettable  case,  the  algorithm  keeps  all  data  samples  and 
the  estimated  frequencies  have  low  variance.  On  the  other  hand,  when  a  is 
approaching  to  0,  which  represents  the  most  forgetful  case,  the  algorithm  keeps 
only  a  few  data  samples  and  the  estimated  frequencies  have  high  variance. 

As  mentioned  above,  the  RTLS  algorithm  applies  the  eigen-decomposition 
technique  to  the  newly  augmented  matrix  C(n).  The  estimated  time-varying 
correlation  matrix  R(n)  is  expressed  as 

R(n)  =  aR(n_D  +  (1  -  a)y(n)y(Hn)  (4.3) 


where  R(n)  is  defined  as 


R(n)  -C(n)C(n). 


(4-4) 


In  this  tracking  (  updating  )  method,  we  use  the  rank-one  eigenstructure 
update  recursively.  In  order  to  diminish  the  computational  load,  we  use  the 
sphericity  property  proposed  by  DeGroat  [12].  Spherical  updating  forces  both 
of  the  eigenvalues  in  the  signal  subspace  and  the  noise  subspace  to  be  replaced 
with  their  mean  values.  Thus,  we  form  two  different  spherical  subspaces 
without  changing  the  dimension  of  the  diagonal  matrix  D  of  R(n).  The  spherical 
subspace  updating  might  be  viewed  as  the  nucleus  of  the  RTLS  adaptive  filter 
algorithm.  We  implement  the  eigendecomposition  of  the  one-step  past-time 
correlation  matrix  R(n-i)  to  get  the  updated  decomposition  at  time  n.  Thus,  we 
have  the  following  updating 

R(n-1)  =  U(n-l)D(n-l)lJ(i_i)  (4.5) 
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(4.6) 


R(„,=U(„)D(„U<H„, 

where  for  all  present  time  n 

«(.)»£, -OS)Uw- 1  (4-7) 

and  D(n)  is  still  a  diagonal  matrix  consisting  of  the  eigenvalues. 

This  updating  algorithm  is  developed  to  track  the  smaller  averaged 
eigenspace  of  the  exponentially  weighted  correlation  matrix  R(n)  under  the 
condition  of  forced  sphericity  of  the  eigenstructure.  To  determine  the  right 
singular  subspace  of  C(n),  we  need  to  obtain  the  eigenspaces  of  R(n) 
corresponding  to  the  right  singular  subspaces  of  C(n).  For  later  simulation,  we 
get  the  vector  y(n)  extracted  from  the  last  row  of  matrix  [  Alb  ].  Next,  for  the 
purpose  of  simplifying  the  derivation  of  the  subspace  updating,  we  neglect  the 
factor  (  1-a  )  and  adopt  a  new  time-varying  correlation  matrix  R  for  the  rank- 
one  updating  recursion 


R(n)  -  otR(n)  +  y(n)y(n) 


By  substituting  (4.6)  into  (4.8)  and  letting 

Z  =  UHyf 


we  further  get 


y  =  ux 


(4.8) 


(4.9) 


(4.10) 
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and  derive  the  following  form 


R  =  U(oD+XXH)IJh  (4.11) 

Since  both  of  the  RTLS  algorithm  and  the  TLS  algorithm  use  the  SVD- 
based  partition  method,  they  have  the  same  solution  for  the  signal  subspace  and 
the  noise  subspace,  respectively.  Namely,  we  can  rewrite  (2.36)  and  (2.39) 
from  Chapter  II  below: 
for  the  signal  subspace,  the  solution  vector  is 


V,g 

xrtls  - ,  mi 
1-g  g 


for  the  noise  subspace,  the  solution  vector  is 


V?c* 

xrtls  ~~~yr 
c  c 


(4.12) 


(4.13) 


In  simulations  considered  later  in  this  chapter,  we  use  two  signal 
frequencies:  one  frequency  is  fixed  and  the  other  one  is  linear  time- varying,  so 
the  smaller  subspace  is  the  signal  subspace.  We  assume  that  these  two 
complex  sinusoidal  signals  have  frequency  components  f2=0.27Hz  and  fx 
moving  around  0.23Hz  with  a  moving  weight  of  3%  and  10%;  the  phase 

7t 

components  are  <)>.=— ,  <(>2=0;  10  independent  trials  and  10  recursive  realizations 

4 

are  used  in  this  simulation.  Note  that  we  have  four  parameters  to  vary:  the 
filter  order  P,  the  exponential  weighting  factor  a,  the  signal-to-noise  ratio 
(  SNR  )  and  the  moving  weight  of  the  time-varying  signal  frequency. 
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3.  Simulations 

First,  we  apply  the  RTLS  algorithm  and  use  the  spherical  subspace 
updating  technique.  Next,  we  compare  the  performance  of  the  RTLS  algorithm 
with  spherical  updating  with  that  of  the  RTLS  algorithm  without  spherical 
updating.  A  data  length  of  100  points  and  a  window  length  of  50  points  are 
used  in  the  following  simulations.  In  addition,  we  assign  the  forgetting  factor  a 
to  be  one  of  the  two  values:  0.25  and  0.95.  Next,  we  adopt  one  of  the  two 
moving  frequency  weights:  3%  and  10%.  We  would  also  like  to  see  when  the 
resolution  of  signal  zeros  breakdowns  by  selecting  various  filter  orders.  Finally, 
by  fixing  a=0.95  and  the  moving  weight=10%,  we  decrease  SNR  to  OdB  and 
-5dB  so  as  to  observe  the  influence  of  low  SNR  on  the  resolution.  The  cases 
studied  are  listed  as  follows  : 

a.  Case  (1) 

Fix  SNR=12dB,  P=8,  moving  weight=3% 

Change  a=0.25, 0.95. 

b.  Case  (2) 

Fix  SNR=12dB,  P=8,  moving  weight=10% 

Use  a=0.95 

c.  Case  (3) 

Fix  SNR=12dB,  a=0.95,  moving  weight=10% 

Choose  P=3,  5. 

d.  Case  (4) 

Fix  P=8,  a=0.95,  moving  weight=10% 

Decrease  SNR=0dB,  -5dB. 
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B.  EXPERIMENTAL  OBSERVATIONS 

1.  Estimated  Signal  Subspace  Obtained  with  Averaging  the 

Eigenvalues 

The  signal  subspace  information  is  used  to  estimate  the  two  signal 
frequencies.  By  using  the  data  in  case  (1),  for  a=0.25,  we  observe  that  the 
locations  of  the  zeros  are  randomly  distributed.  It  is  very  difficult  to  distinguish 
signal  zeros  from  noise  zeros,  as  shown  in  Figure  4.1.1.  The  averaged 
eigenvalues  of  the  spherical  updated  matrix  are  presented  in  Figure  4.1.2.  If  a 
is  chosen  much  higher,  say  a=0.95,  we  observe  that  the  resolution  becomes 
very  good  and  we  can  discern  the  fixed  signal  and  the  moving  signal.  It  is 
noted  that  the  noise  zeros  just  congregate  forming  a  small  cluster  and  are 
distributed  inside  the  unit  circle,  as  shown  in  Figures  4.1.3  and  4.1.4.  Hence,  the 
experiments  show  that  the  case  of  a=0.95  has  much  better  resolution  since  it 
keeps  more  data  samples  as  shown  in  Figure  4.2.1.  Next,  we  observe  that  it  is 
quite  obvious  to  distinguish  the  fixed  frequency  (  a  slight  shift  )  and  the  time- 
varying  frequency  (  countable  small  circles  )  using  the  data  given  in  case  (2). 
The  averaged  eigenvalues  are  also  displayed  in  Figure  4.2.2.  Thus,  results 
show  that  the  case  with  larger  moving  weight  and  higher  forgetting  factor 
(  approaching  to  one  )  has  the  best  zeros  resolution  and  tracking  capability. 
Further  speaking,  the  forgetting  factor  plays  a  more  important  role  than  the 
moving  weight.  By  using  the  data  in  case  (3),  we  observe  that  in  the  case  of 
P=5  the  "  fixed  "  signal  estimate  shows  some  some  variations.  However,  we 
are  still  able  to  differentiate  signal  zeros  from  noise  zeros,  as  shown  in  Figure 

4.3.1.  The  performance  gets  worse  when  we  use  P=3,  as  shown  in  Figure 

4.3.2.  Usually  the  filter  order  number  (  P  )  is  chosen  larger  than  the  true  signal 
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number  ( M  ),  so  we  do  not  consider  the  case  P=M  here.  Using  the  data  in  case 
(4),  we  observe  that  the  resolution  becomes  worse  at  very  low  SNR  since 
some  noise  zeros  are  very  close  to  or  on  the  unit  circle  and  can  be  mistaken  for 
true  signal  zeros,  as  shown  in  Figures  4.4.1  &  4.4.2. 

2.  Estimated  Signal  Subspace  Obtained  without  Averaging 
Eigenvalues 

Next,  we  consider  the  estimated  correlation  matrix  R(n)  obtained 
without  using  averaged  eigenvalues.  Using  the  data  in  case  (2),  we  obtain  very 
good  resolution  of  the  fixed  frequency  and  the  time-varying  frequency,  as 
shown  in  Figure  4.5.1.  The  non-averaged  eigenvalues  of  diagonal  elements  of 
the  correlation  matrix  are  displayed  in  Figure  4.5.2.  Next,  using  the  data  in  case 
(4)  shows  that  the  resolution  becomes  worse  with  very  lower  SNR,  as  shown 
in  Figure  4.6  to  be  compared  with  Figure  4.4.  Besides,  from  the  plot  of  non- 
averaged  eigenvalues,  we  see  the  signal  eigenvalues  are  still  much  higher  than 
the  noise  eigenvalues  as  shown  in  Figures  4.6.2  &  4.6.4. 

3.  Estimated  Noise  Subspace  Obtained  with  Averaging  the 
Eigenvalues 

Here,  we  use  the  spherical  subspace  updating  method  and  choose  the 
noise  subspace  to  estimate  the  frequencies.  The  data  are  generated  using  the 
characteristics  of  case  (2)  given  in  section  3.  The  results  are  quite  different 
from  those  obtained  when  choosing  the  signal  subspace.  Some  of  the  noise 
zeros  appear  on  the  unit  circle  and  can  potentially  be  mistaken  for  signal  zeros, 
as  shown  in  Figure  4.7  to  be  compared  with  Figure  4.2.  Next,  using  the  data 
defined  in  case  (3)  earlier,  we  still  obtain  poor  resolution,  as  shown  in  Figure  4.8 
to  be  compared  with  Figure  4.3.  Next,  we  investigate  the  effects  of  the  SNR  on 
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the  resolution.  Since  the  resolution  of  SNR=12dB  plotted  in  Figure  4.7.1  is  not 
good,  we  raise  the  SNR  to  15dB  and  30dB.  We  observe  that  the  resolution 
does  not  improve,  as  shown  in  Figures  4.9.1  &  4.9.2.  Therefore,  we  conclude 
that  it  is  much  better  and  more  efficient  to  adopt  the  signal  subspace  in  tracking 
two  closely-spaced  signal  sources  (  targets  )  instead  of  adopting  the  noise 
subspace,  when  using  the  spherical  property. 
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V.  CONCLUSIONS 


In  this  thesis,  we  first  introduce  the  Least  Squares  (  LS  )  algorithm  which 
only  considers  the  errors  or  perturbations  in  the  observation  vector  b  of  a  linear 
prediction  equation  Ax=b.  Since  the  assumption  is  not  realistic,  we  present  the 
Total  Least  Squares  (  TLS  )  algorithm  which  considers  the  errors  or 
perturbations  in  the  observation  vector  b  and  the  data  matrix  A.  When  there 
exists  some  linear  dependence  algebraical  relationship  among  the  correlated 
noise  perturbation  components  in  A  and  b,  we  need  to  modify  the  TLS 
algorithm  by  taking  into  the  dependence  account,  which  leads  to  the 
Constrained  Total  Least  Squares  (  CTLS  )  algorithm.  Next,  we  present  the 
Linear  Constraint  Total  Least  Squares  (  LCTLS  )  algorithm.  The  LCTLS  can 
be  viewed  as  a  special  case  of  the  CTLS  algorithm  and  is  used  to  resolve 
closely-spaced  frequencies  with  a  given  known  frequency  component.  Then, 
we  present  numerical  implementations  of  the  TLS  algorithm.  We  investigate 
the  performances  of  the  TLS  algorithm  to  resolve  two  damped  /  undamped 
complex  sinusoidal  signals  in  additive  white  noise.  Next,  we  present  a 
numerical  implementation  of  the  Recursive  Total  Least  Squares  (  RTLS  ) 
algorithm  to  locate  a  fixed  signal  frequency  and  to  track  a  time-varying  signal 
frequency.  The  RTLS  algorithm  uses  an  eigendecomposition  technique  instead 
of  a  singular  value  decomposition  technique.  We  apply  a  spherical  subspace 
updating  method  in  order  to  efficiently  reduce  the  computational  load  and 
compare  the  performances  of  the  RTLS  with  and  without  spherical  updating. 

From  the  simulation,  we  observe  the  case  with  larger  moving  weight  and 
higher  exponential  weighting  factor  (  forgetting  factor  )  has  much  better 
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resolution  and  tracking  capability,  when  the  spherical  signal  subspace  is  used. 
Furthermore,  we  show  that  the  exponential  weighting  factor  plays  a  more 
important  role  than  the  moving  weight.  In  addition,  we  see  that  the  resolution 
becomes  worse  at  very  low  SNR  because  some  noise  zeros  are  approaching  to 
the  unit  circle  and  can  be  potentially  mistaken  for  signal  zeros.  When  we  use 
the  signal  subspace  without  using  averaged  eigenvalues,  we  observe  that  the 
results  do  not  have  any  significant  difference  from  the  case  using  averaged 
eigenvalues.  However,  if  we  use  the  noise  subspace  with  spherical  averaged 
eigenvalues,  poor  resolution  is  obtained.  Therefore,  it  is  better  to  adopt  the 
signal  subspace  for  resolving  two  closely-spaced  signal  sources. 
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TABLE  1.  EFFECTIVE  STANDARD  DEVIATIONS  OBTAINED 
BY  USING  THE  TLS  ALGORITHM 


SNR  (dB) 

Estimated  Deviation  for 

fi=0.23Hz 

Estimated  Deviation  for 

f2=027Hz 

20 

0.0020 

8.5083e-04 

18 

0.0030 

9.1698e-04 

16 

0.0049 

0.0046 

14 

0.0047 

0.0033 

12 

0.0097 

0.0073 

10 

0.0058 

0.0065 

8 

0.0054 

0.0203 

7 

0.0189 

0.0133 

6 

0.0198 

0.0164 

5 

0.0073 

0.0166 

4 

0.0217 

0.0388 

3 

0.0106 

0.0461 

2 

0.0231 

0.0333 

1 

0.0325 

0.0108 

0 

0.0684 

0.1816 
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-1.5 


30  independent  trials 


-1.5 


-1  -0.5  0 
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1  1. 


Figure  3.1.1  TLS  solution  for  SNR=20dB,  P=8 
N=50,  fi=0.23,  f2=0.27,  undamped  sinusoids 


N=50,  fi=0.23,  f2=0.27,  undamped  sinusoids 


Figure  3.1.3  TLS  solution  for  SNR=5dB,  P=8 
N=50,  fi=0.23,  f2=0.27,  undamped  sinusoids 


Figure  3.1.4  TLS  solution  for  SNR=3dB,  P=8, 
N=50,  fi=0.23,  f2=0.27,  undamped  sinusoids 


fl=0.25,  f2=0i27 


N=50,  fi=0.25,  fr=0.27,  undamped  sinusoids 


Figure  3.2.2  TLS  solution  for  SNR=12dB,  P=8 
N=50,  f  1=0.25,  fr=0.27,  undamped  sinusoids 


Figure  3.2.3  TLS  solution  for  SNR=3dB,  P=8.  N=50 
fi=0.25,  f2=0.27,  <|>i=7c/4,  <j>2=-n/4,  undamped  sinusoids 


Figure  3.2.4  TLS  solution  for  SNR=12dB,  P=8,  N=50 
fi=0.25,  f2=0.27,  4>i=tc/4,  4>2=-7t/4,  undamped  sinusoids 


fl=0.25,  f2=0:27 
phil=pi/4,  phi2=0 


10  independent  trials 

Figure  3.3.3  TLS  solution  for  SNR=12dB,  P=8 
N=50,  fi=0.25,  f2=0.27,  undamped  sinusoids 


fl=0.25,  f2=0.27 
phil=pi/4,  phi2=0 


10  independent  trials 


Figure  3.3.4  TLS  solution  for  SNR=12dB,  P=18 
N=50,  fi=0.25,  f2=0.27,  undamped  sinusoids 


1.5 


Figure  3.4.1  TLS  solution  for  SNR=12dB,  P=8 
N=50,  fi=0.25,  f2=0.27,  undamped  sinusoids 
signal  and  noise  subspaces  mixed  by  3  columns 


Figure  3.4.2  TLS  solution  for  SNR=12dB,  P=8 
N=50,  fi=0.25,  f2=0.27,  undamped  sinusoids 
signal  and  noise  subspaces  mixed  by  6  columns 
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N=50,  fi=0.25,  £2=0.27,  damped  sinusoids 


N=50,  fi=0.25,  £2=0.27,  damped  sinusoids 
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alphal=0.8 

alpha2=0.9 


fl=0.25,  f2=0.27 
phil=pi/4,  phi2=0 


30  independent  trials 

Figure  3.5.3  TLS  solution  for  SNR=5dB,  P=18 
N=50,  fi=0.25,  ft=0.27,  damped  sinusoids 


alphal=0.8 

alpha2=0.9 


fl=0.25,f2=0.27 
phil=pi/4,  phi2=0 


[HE 


;  30  independent  trials 
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-1.5  -1  -0.5  0  0.5  1  1.5 

Figure  3.5.4  TLS  solution  for  SNR=20dB,  P=18 
N=50,  fi=0.25,  f2=0.27,  damped  sinusoids 
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Figure  4.1.1  RTLS  solution  for  SNR=12dB,P=8,N=100 
<x=0.25,  weight=3%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’S’ 


Figure  4.1.2  averaged  eigenvalues  of  updated  spherical  matrix  'S' 
SNR=12dBJ>=8,N=100,a=0.25,weight=3%,fi  is  time-varying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.1.3  RTLS  solution  for  SNR=12dBrP=8,N=100 
Ot=0.95,  weight=3%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’S' 
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Figure  4.1.4  averaged  eigenvalues  of  updated  spherical  matrix  'S’ 
SNR=12dB  J>=8,N=100,a=0.95,weight=3%,fi  is  time-varying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.2.1  RTLS  solution  for  SNR=12dB,P=8,N=100 
a=0.95,  weight=10%,  fl  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’S’ 


Figure  4.2.2  averaged  eigenvalues  of  updated  spherical  matrix  ’S' 
SNR=12dB,P=8,N=100,a=0.95,weight=10%,fi  is  time-varying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.3.1  RTLS  solution  for  SNR=12dB,  P=5,  N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’S’ 


Figure  4.3.2  RTLS  solution  for  SNR=12dB,  P=3,  N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S’ 
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Figure  4.4.1  RTLS  solution  for  SNR=OdBJ>=8,N=100 
a=0.95,  weigh t=  10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ‘S' 
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Figure  4.4.2  RTLS  solution  for  SNR=-5dB,P=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’S’ 
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Figure  4.5.1  RTLS  solution  for  SNR=12dB,P=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’R’ 
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Figure  4.5.2  non-averaged  eigenvalues  of  correlation  matrix  'R' 
SNR=12dB,P=8,N=100,a=0.95,weight=10%,fi  is  time-varying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.6.1  RTLS  solution  for  SNR=OdBJP=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’R’ 
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Figure  4.6.2  non-averaged  eigenvalues  of  correlation  matrix  'R' 
SNR=0dB,P=8,N=100,a=0.95,weight=10%,fi  is  time-vaiying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.6.3  RTLS  solution  for  S  NR=-5dB  ,P=8 ,N=  1 00 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  signal  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  ’R* 


Figure  4.6.4  non-averaged  eigenvalues  of  correlation  matrix  'R' 
SNR=-5dB,P=8,N=100,Ot=0.95,weight=10%,fl  is  time-varying 
solution  derived  by  using  the  signal  subspace  information 
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Figure  4.7.1  RTLS  solution  for  SNR=12dB,P=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  noise  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S' 


Figure  4.7.2  averaged  eigenvalues  of  updated  spherical  matrix  ’S’ 
SNR=12dB  J>=8,N=100,(X=0.95,weight=10%,fi  is  time-varying 
solution  derived  by  using  the  noise  subspace  information 
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Figure  4.8. 1  RTLS  solution  for  SNR= 1 2dB  J*=5,N=1 00 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  noise  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S' 


Figure  4.8.2  RTLS  solution  for  SNR=12dB,P=3,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  noise  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S’ 
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Figure  4.9.1  RTLS  solution  for  SNR=15dBJP=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  noise  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S' 


Figure  4.9.2  RTLS  solution  for  SNR=30dB,P=8,N=100 
a=0.95,  weight=10%,  fi  is  time-varying,  using 
the  noise  subspace  information  &  performing 
the  eigendecomposition  of  spherical  matrix  'S’ 
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APPENDIX  A  :  MATLAB  SOURCE  CODE 


This  appendix  consists  of  MATLAB  codes  used  to  implement  the  Total 
Least  Squares  (  TLS  )  algorithm  mentioned  in  Chapter  III  and  the  Recursive 
Total  Least  Squares  (  RTLS  )  algorithm  mentioned  in  Chapter  IV. 

1.  IMPLEMENTATION  OF  THE  TLS  ALGORITHM 

clear 

%initialize  the  set  of  roots 
rt=[]; 

nn=input(’  trial  number '); 

SNR=input(’  S/N(dB)  =  '); 
p=input('  prediction  filter  order  '); 
y  1  =zeros(p,  1  );y2=zeros(p,  1 ); 

%construct  the  received  signal  and  begin  the  trial 
for  h=l  :nn. 


n=0:49; 

randCnormal'); 

wl=rand(l,50); 

randCnormal'); 

w2=rand(l,50); 

f  1=0.25; 

%signal  frequencies 

f2=0.27; 

%f  1=0.23  for  option 

phil=pi/4; 

%signal  phase  angles 

%decaying  constants 


alphal=0.8; 

alpha2=0.9; 

j=sqrt(-l); 

wn=wl+j*w2; 

sl=exp(j*(2*pi*fl  *n+phil)); 
s2=exp(j  *(2*pi*f2*n+phi2)); 
s  l=exp(-alphal  *n+j*(2*pi*f  1  *n+phil )); 
s2=exp(-alpha2*n+j*(2*pi*f2*n+phi2)); 
ss=sl+s2; 

xn=ss+sqrt(l/(2*10A(SNR/10))).*wn; 
nsin=2; 
n=50; 
if  n>p+2, 

for  k=l  :n-p, 

A(k,  1  :p)=xn(k:p+k- 1 ); 
end 

b=xn(l,p+l:n); 

bb=conj(b)'; 

end 

C=[bb  A]; 

[u,s,v]=svd(C); 

sigma=diag(s); 

ct=v(l,nsin+l:p+l); 

c=conj(ct)'; 

vp= v(2:p+ 1  ,nsin+ 1  :p+ 1 ); 


%independent  complex  noise 
%create  two  undamped  signals 

%create  two  damped  signals 


%received  resulted  signal 
%for  two  sinusoidal  sources 
%window  data  sample  number 


%generate  data  matrix 

%generate  observation  vector 


%create  augmented  matrix 
%perform  singular  value  decom. 

%extract  the  desired  vector 
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v2=[ct;vp]; 

xtls=-(vp*conj(c))./(conj(ct)*c);  %solution  from  noise  subspace 

z=[l  -fliplr(conj(xtls)’)];  %prediction  filter  polynomial 

rt=[rt  roots(z)]; 

zl=  0.3+j*0.7;  %create  sector  regions  vs.  origin 

z2=-0.3+j*0.7; 

real_rt=real(rt); 

for  nl=l:p,  %calculate  estimated  deviations 

if  real_rt(nl,h)  >  0, 
yl(nl,h)=rt(nl,h); 
yll=sort(yl); 
yll=flipud(yll); 
new_y  11(1  :p/2,h)=y  1 1  ( 1  :p/2,h); 
kl=find((angle(new_yll(l:p/2,h))  >  angle(zl))... 

&  (angle(new_y  11(1  :p/2,h)>  <  pi/2»; 
if  length(kl)~=l, 

error('More  than  1  element  of  new_yll  is  nonzero.’); 
end 

esti__yl  l(l,h)=new_yl  l(kl,h); 
estdl=sqrt(l/nn*(sum(abs(esti_y  1  l(l,h)... 

-(0.1253+j *0.9921)  .*ones(l,nn))  .A2))); 

else 

y2(nl,h)=rt(nl,h); 

y22=sort(y2); 

y22=flipud(y22); 
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new_y22(  1  :p/2,h)=y22(  1  :p/2,h); 
k2=find((angle(new_y22(l:p/2,h))  >  pi/2)... 

&  (angle(new_y22(l:p/2,h))  <  angle(z2))); 
if  length(k2)~=l, 

error('More  than  1  element  of  new_y22  is  nonzero.'); 
end 

esti__y22(l  ,h)=new_y22(k2,h); 
estd2=sqrt(l/nn*(sum(abs(esti_y22(  l,h)... 

-(-0.1253+j*0.9921)  .*ones(l,nn))  .A2))); 
end 
end 
end 

%plot  the  zeros  distribution  situations 

axis(’square') 

axis([-1.5  1.5 -1.5  1.5]) 

plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),’-w',... 

real(new_y  1 1  ),imag(new_y  1 1  ),,o\real(new_y22),imag(new_y22),,o'); 

grid 

title('  TLS  solution  for  SNR=  ,  P=  ,  N=  ’); 

text(-l.l,1.3,'  alphal=  '); 

text(- 1 . 1 , 1 . 1 alpha2=  '); 

text(0,1.3,'  fl=  ,  f2=  '); 

text  (0,1.1,'  phil=  ,  phi2=  '); 

text(0,-1.3,'  30  independent  trials  '); 

text(-l,-1.4,'  destroyed  signal  subspace  by  #  columns  ’ ); 
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%REMARKS: 

%(l).When  SNR  down  to  1  dB,  we  need  change  z2  to  -0.4+j*0.8 
%  in  order  to  make  k2  for  new_y22  satisfy  'find'  condition. 
%(2).For  SNR  down  to  0  dB,  we  have  to  change  zl  to  0.4+j*0.7 
%  in  order  to  make  kl  for  new_yl  1  satisfy  'find'  condition. 
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2.  IMPLEMENTATION  OF  THE  RTLS  ALGORITHM  USING 
THE  SIGNAL  SUBSPACE  AND  PERFORMING  THE  EIGEN- 
DECOMPOSITION  OF  MATRIX  '  S  ’ 

Note  that  S  =  a*D+(l-tx)*3*3H  3  =  UH*x 

clg 

clear 

%initialization  of  conditions 
rt=[];esp=[];sig=[];noi=[]; 
nt=input('  trial  number= '); 
ni=input('  iteration  number= '); 
p  =input('  prediction  order=  '); 

SNR=input('  signal-to-noise  ratio,  S/N(dB)=  '); 
weight=input('  the  weight  of  moving  frequency=  '); 
alpha=input(’  exponential  weighting  factor,  alpha=  ’); 

%construct  received  signal  and  begin  independent  trial 
for  h=l:nt, 
nl=l:100; 
rand('normal’); 
wl=rand(l,100); 
rand('normal'); 
w2=rand(l,100); 

f  1  =0.23-(h- 1  )*weight/nt;  %make  fi  time-varying  moving 

f2=0.27; 

phi  1  =pi/4; 

phi2=0; 
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^independent  complex  noise 
%create  two  undamped  signals 


j=sqrt(-l); 

wn=wl+j*w2; 

s  1  =exp(j  *(2*pi*f  1  *nl  +phi  1 )); 
s2=exp(j  *(2*pi*f2*nl  +phi2)); 
ss=sl+s2; 

xn=ss+sqrt(l/(2*10A(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if  n2>p+2, 

for  k=l  :n2-p, 

A(k,  1  :p)=xn(k:p+k- 1 ); 
end 

b=xn(l,p+l:n2); 

bb=conj(b)'; 

end 

B=[bb  A]; 
aph=alpha  .Anl; 
aph=fliplr(aph(l  :n2-p)); 
T=sqrt(diag(aph)); 

C=T*B; 

R=C'*C; 

[U,D]=eig(R); 

fort=l:ni. 


%received  resulted  signal 
%for  two  sinusoidal  sources 
%  window  data  sample  number 


%generate  data  matrix 


%generate  observation  vector 


%create  augmented  matrix 
%define  estd.  correlation  matrix 
%perform  eigendecomposition 
%recursively  update  subspace 


y=xn(  1  ,n2-p+t:n2+t); 


beta=U'*x; 

D 1  =sort((diag(D))); 

D2=flipud(Dl); 

D2=D2'; 

s_eigval=mean(D2(l  :nsin)); 
n_eigval=mean(D2(nsin+ 1  :p+l)); 
s_eigspa=s_eigval  .*eye(nsin); 
n_eigspa=n_eigval  .*eye(p-nsin+l); 
s_diaspa=diag(s_eigspa); 
n_diaspa=diag(n_eigspa); 

Dhat  1 =[s_eigspa  zeros(nsin,p-nsin+ 1 )] ; 

Dhat2=[zeros(p-nsin+l,nsin)  n_eigspa]; 

Dhat  =[Dhatl;Dhat2]; 

S=alpha  .*Dhat+(l -alpha)  .*beta*beta';  %our  major  concern 

[UU  ,DD]=eig(S); 

Unew=U*UU; 

Dnew=DD; 

gt=Unew(l,l:nsin);  %extract  the  desired  vector 

g=conj(gt)’; 

u  lp=Unew(2:p+ 1 , 1  :nsin); 
ul=[gt;ulp]; 

xtls=(ulp*conj(g))./(l-conj(gt)*g);  %solution  from  signal  subspace 

z=[l  -fliplr(conj(xtls)’)]; 
rt=[rt  roots(z)]; 
esp=[esp  D2']; 
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sig=[sig  s_diaspa]; 
noi=[noi  n_diaspa]; 

U=Unew;  %assign  new  U  &  D  for  updating 

D=Dnew; 
end 
end 

%plot  zeros  in  complex  z-plane 

axis('square'); 

axis([-1.5  1.5  -1.5  1.5]); 

plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),'-w',... 

real(rt),imag(rt),'o'),grid 

title(’  RTLS  solution  for  SNR=  dB,  P=  ,  N=  ’); 

text(-l,1.3,'  alpha=  '); 

text(- 1 ,1.1,'  weight=  '); 

text(0. 1,1.1,'  fl  moving,  f2=0.27  '); 

text(0.1,1.3,'  phil=pi/4,  phi2=0 '); 

text(-0.5,-1.2,'  10  independent  trials  '); 

text(-0.5,-1.4,’  10  recursive  realizations  '); 

%plot  spherical  singularvalues  status  of  updated  matrix 

clg 

axis; 

n3=l:p+l; 
plot(n3  ,esp,'o'),grid 
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3.  IMPLEMENTATION  OF  THE  RTLS  ALGORITHM  USING 
THE  SIGNAL  SUBSPACE  AND  PERFORMING  THE  EIGEN- 
DECOMPOSITION  OF  MATRIX  ’  R  ’ 

Note  that  R  =  a*R+(l-a)*x*xH 
clear 

%inidalization  of  conditions 
rt=[];esp=[]; 

nt=input(’  trial  numbers '); 
ni=input('  iteration  numbers  ’); 
p  =input(’  prediction  orders  '); 

SNR=input('  signal-to-noise  ratio,  S/N(dB)=  '); 
weight=input('  the  weight  of  moving  frequency=  '); 
alpha=input('  exponential  weighting  factor,  alphas  '); 

%construct  received  signal  and  begin  independent  trial 
for  h=l:nt, 
nl=l:100; 
rand('normal'); 
wlsrand(  1,100); 
randCnormal'); 
w2=rand(  1,100); 

f  1  =0.23-(h- 1  )*weight/nt;  %make  f l  time- varying  moving 

f2=0.27; 

philsf  (4; 

phi2=0; 

j=sqrt(-l); 
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wn=wl+j*w2; 

s  1  =exp(j  *(2*pi*f  1  *n  1  +phi  1 )); 
s2=exp(j  *(2*pi*f2*n  1  +phi2)); 
ss=sl+s2; 

xn=ss+sqrt(l/(2*10A(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if  n2>p+2, 

for  k=l  :n2-p, 

A(k,  1  :p)=xn(k:p+k- 1 ); 
end 

b=xn(l,p+l:n2); 

bb=conj(b)’; 

end 

B=[bb  A]; 
aph=alpha  .Anl; 
aph=fliplr(aph(l  :n2-p)); 
T=sqrt(diag(aph)); 

C=T*B; 

R=C’*C; 

[U,D]=eig(R); 

fort=l:ni, 

y=xn(  1  ,n2-p+t:n2+t); 
x=y'; 

R=alpha  .*R+(1 -alpha)  .*x*x'; 


%independent  complex  noise 
%create  two  undamped  signals 

%received  resulted  signal 
%for  two  sinusoidal  sources 
%window  data  sample  number 

%generate  data  matrix 
%generate  observation  vector 


%create  augmented  matrix 
%define  estd.  correlation  matrix 
%perform  eigendecomposition 
%recursively  update  subspace 

%our  major  concern 
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[UU,DD]=eig(R); 

D 1  =flipud(sort((diag(DD))»; 

gt=UU(l,l:nsin);  ^extract  the  desired  vector 

g=conj(gt)'; 

ulp=UU(2:p+l,l:nsin); 

ul=[gt;ulp]; 

xtls=(ulp*conj(g))./(l-conj(gt)*g);  %solution  from  signal  subspace 

z=[l  -fliplr(conj(xtls)')]; 
rt=[rt  roots(z)]; 
esp=[esp  Dl]; 
end 
end 

%plot  zeros  in  complex  z-plane 

axis(’square'); 

axis([-1.5  1.5  -1.5  1.5]); 

plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),’-w\... 

real(rt),imag(rt),,o'),grid 

titleC  RTLS  solution  for  SNR=  dB,  P=  ,  N=  ’); 
text(-l,1.3,’ alpha=  ’); 
text(- 1,1.1,'  weight= '); 
text(0. 1,1.1/  fl  moving,  f2=0.27  '); 
text(0.1,1.3,’ phil=pi/4,  phi2=0  ’); 
text(-0.5,-1.2,'  10  independent  trials  '); 
text(-0.5,-1.4,'  10  recursive  realizations  ’); 
clg 
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axis; 

n3=l:p+l; 

plot(n3,esp,'o'),grid 


4.  IMPLEMENTATION  OF  THE  RTLS  ALGORITHM  USING 
THE  NOISE  SUBSPACE  AND  PERFORMING  THE  EIGEN- 


DECOMPOSITION  OF  MATRIX  1  S  ’ 

Note  that  S  =  a*D+(l-a)*p*J3H  (5  =  UH*x 

clg 
clear 

%inidalization  of  conditions 
rt=[];esp=[];sig=[];noi=[]; 
nt=input('  trial  number= '); 
ni=input(’  iteration  number= '); 
p  =input(’  prediction  order=  '); 

SNR=input('  signal-to-noise  ratio,  S/N(dB)=  ’); 
weight=input('  the  weight  of  moving  frequency=  ’); 
alpha=input(’  exponential  weighting  factor,  alpha=  '); 

%construct  received  signal  and  begin  independent  trial 
forh=l:nt, 
nl=l:100; 
rand('normal'); 
wl=rand(l,100); 
randCnormal’); 
w2=rand(l,100); 

f  1  =0.23-(h- 1  )*weight/nt;  %make  fi  time-varying  moving 

f2=0.27; 

phil=pi/4; 

phi2=0; 
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^independent  complex  noise 
%create  two  undamped  signals 


j=sqrt(-l); 

wn=wl+j*w2; 

s  1  =exp(j  *(2*pi*f  1  *n  1  +phi  1 )); 
s2=exp(j  *(2*pi*f2*n  1  +phi2)); 
ss=sl+s2; 

xn=ss+sqrt(l/(2*10A(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if  n2>p+2, 

for  k=l  :n2-p, 

A(k,  1  :p)=xn(k:p+k- 1 ); 
end 

b=xn(l,p+l:n2); 

bb=conj(b)'; 

end 

B=[bb  A]; 
aph=alpha  .Anl; 
aph=fliplr(aph(  1  :n2-p»; 
T=sqrt(diag(aph)); 

C=T*B; 

R=C’*C; 

[U,D]=eig(R); 

fort=l:ni, 

y=xn(  1  ,n2-p+t:n2+t); 
x=y'; 


%received  resulted  signal 
%for  two  sinusoidal  sources 
%window  data  sample  number 


%generate  data  matrix 


%generate  observation  vector 


%create  augmented  matrix 
%define  estd.  correlation  matrix 
%perform  eigendecomposition 
%recursively  update  subspace 
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beta=lT*x; 

D 1  =sort((diag(D))); 

D2=flipud(Dl); 

D2=D2'; 

s_eigval=msan(D2(  1  :nsin)); 
n_eigval=mean(D2(nsin+l  :p+l)); 
s_eigspa=s_eigval  .*eye(nsin); 
n_eigspa=n_eigval  .*eye(p-nsin+l); 
s_diaspa=diag(s_eigspa); 
n_diaspa=diag(n_eigspa); 

Dhat  1 = [s_eigspa  zeros(nsin,p-nsin+ 1 )] ; 

Dhat2=  [zeros(p-nsin+ 1  ,nsin)  n_eigspa] ; 

Dhat  =[Dhatl;Dhat2]; 

S=alpha  .*Dhat+(l -alpha)  .*beta*beta‘;  %our  major  concern 
[UUJDD]=eig(S); 

Unew=U*UU; 

Dnew=DD; 

ct=Unew(l,nsinH-l:p+l);  %extract  the  desired  vector 

c=conj(ct)'; 

vp=Unew(2:p+ 1  ,nsin+ 1  :p+ 1 ); 
v2=[ct;vp]; 

xtls=-(vp*conj(c))./(l-conj(ct)*c);  %solution  from  noise  subspace 

z=[l  -fliplr(conj(xtls)’)]; 
rt=[rt  roots(z)]; 
esp=[esp  D2’]; 
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sig=[sig  s_diaspa]; 
noi=[noi  n_diaspa]; 

U=Unew;  %assign  new  U  &  D  for  updating 

D=Dnew; 
end 
end 

%plot  zeros  in  complex  z-plane 

axis('square'); 

axis([-1.5  1.5  -1.5  1.5]); 

plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),’-w',... 

real(rt),imag(rt),,o’),grid 

titleC  RTLS  solution  for  SNR=  dB,  P=  ,  N=  '); 
text(- 1,1.3,'  alpha=  '); 
text(- 1,1.1,’  weight= '); 
text(0.1,l.l,'  fl  moving,  f2=0.27  ’); 
text(0.1,1.3,' phil=pi/4,  phi2=0'); 
text(-0.5,-1.2,'  10  independent  trials  '); 
text(-0.5,-1.4,'  10  recursive  realizations  '); 

%plot  spherical  singularvalues  status  of  updated  matrix 

clg 

axis; 

n3=l:p+l; 
plot(n3  ,esp,’o’),grid 
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